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ABSTRACT 

The  development  of  two-  and  three-dimensional  Euler  correction  methods  based 
on  the  Clebsch  transformation  Is  described.  In  these  methods,  the  velocity 
field  Is  decomposed  Into  Irrotatlonal  and  rotational  parts.  A  multigrid  full- 
potential  method  based  on  both  the  finite-difference  and  f Inljte-vclume  formu¬ 
lations  Is  modified  to  solve  for  the  Irrotatlonal  part,  while  the  momentum 
equation  Is  applied  to  solve  for  the  rotational  part.  Two  approaches  are 
developed  to  solve  for  the  rotational  field.  The  approximate  Euler-Clebsch 
approach  assumes  the  entropy  Is  convected  along  mesh  lines,  while  the  exact 
Euler-Clebsch  approach  solves  for  the  convection  of  entropy  along  streamlines. 
The  two  approaches  agree  well  In  the  airfoil  application.  Only  the  approximate 
Euler-Clebsch  approach  Is  employed  In  the  three-dimensional  calculations.  A 
studj  of  finite-difference  and  finite-volume  formulations  of  the  full-potential 
equation  Is  also  Included.  Euler-correction  solutions  are  presented  for  vari¬ 
ous  airfoils,  wings  and  an  F-14  wing/body  and  are  compared  with  results  of  the 
full-potential  and  the  time-marching  Euler  methods. 


1.0  INTROOUCTION 

Invlsdd  flow  generally  can  be  described  by  the  Euler  equations.  However,  the 
use  of  the  Euler  equations  requires  the  flow  density,  velocity  and  energy  to 
be  solved  as  unknown  functions.  Transonic  flows  with  relatively  weak  stacks 
can  be  assumed  to  be  Isent  :p1c  and  approximated  by  potential  flow  methods 
where  the  potential  function  Is  solved  as  the  only  unknown  function,  thus 
requiring  much  less  computation  than  the  Euler  methoas.  In  the  past  fifteen 
years  transonic  computational  methods  for  solving  for  the  potential  flow- 
fields  have  been  well  developed,  extensively  validated,  and  widely  accepted  as 
a  routine  tool  for  aerodynamic  design.  Methods  for  solving  the  Euler  equations 
were  Introduced  at  a  later  stage,  and  significant  progress  has  been  achieved 
In  this  area  recently.  However,  because  of  the  much  greater  memory  and  compu¬ 
tational  time  requirements  of  the  Euler  methods,  and  also  their  less  complete 
validation  [1-2],  the  f ul 1 -potentia 1  methods  are  still  preferred  In  the  routine 
transonic  aerodynamic  design. 

Although  the  existing  full-potential  methods  are  limited  to  problems  dealing 
with  relatively  weak  shocks,  they  are  Ideal  for  simulating  transonic  flowflelds 
0043h  1 


at  cruise  conditions.  The  solution  accuracy  of  the  full-potential  method 

deteriorates  as  the  flow  condition  approaches  stall  where  the  shocks  become 

much  stronger.  To  further  Improve  existing  ful 1 -potential  methods,  and  extend 

their  applications  to  stronger  shock  cases,  the  effects  of  the  total  pressure 

loss  across  shocks  and  the  vortldtles  generated  downstream  of  the  shocks  on 

the  solutions  have  to  be  studied  carefully.  This  can  be  accomplished  by  the 

»  • 

so-called  Euler  correction  methods  that  are  derived  from  the  basic  full- 
potential  methods  by  add-  1ng  additional  terms  to  model  the  nonisentropic 
effects . 

This  report  describes  several  Improvements  of  the  existing  transonic  full- 
potential  method.  Including  the  development  of  an  Euler  correction  method 
using  the  Clebsch  transformation.  Both  finite-difference  and  f Inite-volume 
formulations  are  used  to  solve  the  full-potential  equation  for  the  Irrotatlonal 
part  of  the  flowfleld.  The  finite-difference  formulation  Is  based  on  a  gen¬ 
eralized  coordinate  transformation  [3,4],  while  the  finite-volume  formulation 
Is  based  on  mass  flux  balance  [5].  General  nonconservative,  partially  con¬ 
servative  and  fully  conse-vatlve  artificial  viscosities  and  shock-point  oper¬ 
ators,  as  described  in  Ref.  3,  are  applied  to  reflect  the  directional  bias  of 
local  supersonic  flows.  The  rotational  part  of  the  flowfleld  Is  determled  from 
the  momentum  equation. 

brief  ove'vlew  o:  Euler  correction  methocs  is  given  In  Section  2.0.  The 
f 1n1 te-dl f ference  and  finite-volume  formulations  of  the  full-potential  equa¬ 
tion  are  described  In  Section  3.0’.  The  Clebsch  formulation  of  the  ro*ational 
velocity  components  Is  described  1n  Section  4.0.  The  nume.lcal  procedure  and 
results  for  two-  and  three-dimensional  flows  are  described  In  Sections  5.0  and 
6.0,  respectively.  Concluding  remarks  are  given  In  Section  7.0. 
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2.0  OVERVIEW 


The  Euler  correction  methods  which  produce  Euler-like  solutions  can  be  classi¬ 
fied  Into  three  categories  according  to  their  degree  of  approximation.  The 
first  simply  accounts  for  the  total  pressure  and  density  losses  across  the 

shocks.  The  total  pressure  Is  computed  downstream  of  the  shock  and  remains 

*  * 

constant  along  streamlines.  Because  of  a  total  density  Jump  across  the  shock, 
the  governing  equation  at  shock  points  Is  significantly  modified,  while  the 
governing  equation  downstream  of  the  shock  point  Is  slightly  modified  due  to 
the  variation  of  total  density  between  streamlines.  The  total  pressure  down¬ 
stream  of  the  shock  Is  computed  according  to  the  Rankine-Hugonlot  shock  rela¬ 
tion,  and  no  additional  equation  is  needed  for  Its  solution.  Because  the  vor- 
tlclty  downstream  of  shocks  has  only  third-order  effects  on  the  solutions  In 
the  transonic  region,  this  approach  generally  predicts  shocks  agree'ng  well 
with  Euler  solutions  than  either  full  conservative  or  nonconservative  schemes. 

Several  Investigators  have  developed  schemes  In  this  category,  including  Hafez 
and  Lovell  [6],  and  Kloofer  and  Nixon  [7}.  The  partially  conservatWe  schemes 
developed  by  Lock  [8]  and  Chen  [9]  are  also  Included  In  this  category,  because 
the  addition  of  a  mass  source  at  shock  points  In  the  partially  conservative 
schemes  is  similar  to  the  correction  of  total  density  downstream  of  the  shock. 
In  addition,  the  tota'  pressure  loss  downstream  of  the  shock  can  significantly 
affect  the  Kutta  conqitl'rn.  For  example,  Chen,  C'ark  and  vassberg  [10]  showed 
that  by  imposing  the  noni sentropic  Kutta  condition,  the  comruted  shocks  agree 
fairly  well  with  the  Euler  solutions,  even  for  strong  shock  cases. 

The  second  kind  of  Euler  correction  method  decomposes  the  velocity  vector  Into 
potential  and  rotational  components.  The  rotational  component  Is  explicitly 
related  to  the  vorticity  field  downstream  of  the  shock.  There  are  different 
ways  to  decompose  the  velocity  vector.  Brown.  Brecht  and  Walsh  [11]  simply 
aid  a  scalar  function  to  the  streatmvise  velocity  component.  Sokhey  [12]  and 
hafez  and  Lovell  [6]  apply  Helmholtz  theory  to  decompose  the  velocity  vector 
Into  Irrotatlonal  and  the  rotational  components  such  that  the  stream  function 
can  be  used  to  compute  tne  vorticity.  Ecer  and  Akay  [13,14]  apply  the  Clebsch 
transformation  to  define  the  rotational  velocity  component  In  conjunction  with 
the  use  of  a  finite-element  method.  This  kind  of  Euler  correction  method  Is 
more  general  than  the  first  method  and  additional  governing  equations  are 
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applied  to  solve  for  the  rotational  velocity  components.  The  additional  equa¬ 
tions  for  the  rotational  part  which  govern  the  convection  of  the  vortldty  are 
hyperbolic,  and  require  a  relatively  small  amount  of  computational  time  to 
solve. 

The  third  kind  of  Euler  correction  method,  proposed  by  T.C.  Ta1  [15],  Is  a 
hybrid  method  which  combines  a  finite-difference  relaxation  method  with  the 
method  of  integral  relations.  The  Euler  equations  In  Integral  form  are  solved 
downstream  of  the  shock.  Both  the  shock  location  and  circulation  are  continu¬ 
ously  updated  during  the  Iterations  until  the  far  downstream  and  the  Kutta 
condition  are  satisfied.  Most  of  the  works  cited  above  are  for  two-dimensional 
analysis,  except  the  works  of  Ecer  [14], 

The  present  method  applies  the  Clebsch  transformation  as  the  method  of  Ecer  and 
Akay,  except  that  1r.  the  present  method  the  f Inite-dlf ference/f Inite-volume 
fu  I -potential  method  Is  applied  to  solve  for  the  Irrotatlonal  part  and  also 
that  different  Clebsch  variables  are  chosen.  Both  the  two-  and  three- 
dimensional  full-potential  methods  are  well  developed,  and  the  extension  of 
these  methods  to  include  the  rotational  velocity  components  Is  straightforward, 
furthermore,  the  multigrid  scheme  developed  for  the  full-potential  methods 
works  equally  wel'  in  the  present  method  despite  the  additional  equations  for 
’•oiat^onal  velocity  compon.‘nts. 

A  preliminary  description  of  the  present  me  hod  and  results  is  given  in  Ref. 

16.  Both  the  two-  and  three-dimensional  methods  are  described  In  more  detail 
In  the  following  sections.  Solutions  are  presented  for  various  airfoils,  wings 
and  winq/bodles,  Including  an  F-14  wing/body,  and  comparisons  are  made  with 
both  the  full  potential  and  time  marching  Euler  solutions  when  available. 
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3.0  FINITE-VOLUME  VERSUS  FINITE-DIFFERENCE  FORMULATION 


Both  the  finite-difference  formulation  of  Chen  [3,4]  and  the  finite-volume 
formulation  of  Jameson  and  Caughey  [5]  are  used  In  the  present  method  to  solve 
the  full-potential  equation  for  the  Irrotatlonal  part  of  the  flowfleld.  In 
the  finite-difference  formulation,  the  full-potential  equation  written  In  gen¬ 
eral  coordinates  Is  employed;  since  a  quadratic  coordinate  transformation  Is 
applied,  the  odd-  and  even-point  solutions  are  naturally  coupled  together.  In 
the  finite  volume  formulation,  the  density  term  Is  explicitly  computed  from 
the  potential  function  distribution,  and  fluxes  through  cell  surfaces  are  con¬ 
served.  However,  coupling  terms  are  needed  to  add  to  the  finite  volume  formu¬ 
lation  to  avoid  the  decoupling  of  odd-  and  even-point  solutions.  The  finite- 
difference  formulation  Is  described  1n  Section  3.1,  the  finite-volume  formula¬ 
tion  In  Section  3.2,  and  artificial  viscosities  and  shock-point  operators  are 
descrloed  In  Section  3.3. 

3 . 1  f Inite-Olfference  Method 

The  full-potential  equation  can  be  expressed  as 

(a^  -  >  (a^  -  ♦  (a^  -  "  ^“''^xy  ‘ 

where  u,  v,  w  are  the  x,  y,  z  components  of  the  flow  velocity,  respectively, 
and  a  Is  the  local  speed  of  sound  determined  from  the  energy  equation 

a  *  -  -‘-j—  (u  ♦  V  w  )  (2) 

where  y  Is  the  ratio  of  specific  heats  for  the  assumed  calorically  perfect 
gas  an  a^  Is  the  stagnation  speed  of  sound. 

After  performing  the  matrix  Inversion,  multiplication,  and  algebraic  manipula- 
tion,  a  transformed  full-potential  equation  multiplied  by  the  determinant  of 
the  Jacobian  transformation  matrix,  0,  In  general  curvilinear  coordinates  can 
be  derived  as  [3,4] 

CO  ♦■CO  ♦CO  *00  *00  +00  ♦cO  ♦cO  ♦  c.o  =  0  (3) 

rxx  2^YY  3  ZZ  4^XY  5  YZ  6^XZ  7^X  8^Y  9  2  ^  ' 

where 

2  2  2  2  2 

c^  3  [a  (h^  ♦  h^  ♦  h^)  -  U  ]/0  (4) 
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c  k 
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h;I* 

.*♦5 

:5s 

:* 


Cj  -  [a^(hj  .  ♦  h^)  -  V^]/0 

2  2  2  2  2 
C3  .  [aNh‘  ♦  h‘  ♦  h‘)  -  u^]/0 


>  [2a  (h^h^  ♦  hjhj  ♦  h^h^)  -  2UV]/D 


C5  =  [2a  (h^h^  ♦  hjhg  ♦  h^hg)  -  2VWJ/0 


Cg  =.  [2a  (h^h^  ♦  h^hg  ♦  h^h^)  -  2UW]/0 

c,  .  -  (h,p^  .  .  hjP^j/D 

'8  ■  -  *  \’‘V  ''^ 

S  ■  -  •S*’*  •  Vy  •  N'2>^“ 

L),  V,  W  are  velocity  components  defined  as 

U  =  h,u  ♦  h,v  ♦  h.w 
1  2  3 

V  a  h  u  ♦  h  V  ♦  h  w 
4.  0 

W  a  h,u  ♦  hgV  *  h^w 

Coefficients  h-]  . . .  hg  are  transformation  derivatives  oefined  as 


^^z  - 

'zS 

(■6) 

- 

'Z*Y 

-  ) 

^Y^Z  • 

*Z^ 

(18' 

'z^  - 

(19) 

^Z^X  ■ 

^X^Z 

(20) 

*2^  - 

^x^z 

(21) 

^xS 

^^x 

(22) 

- 

(23) 

^X^Y 

*Y^ 

(24) 
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and  coefficients  Pj^.  Py,  and  p^  are  second-order  transfornatlon 
derivatives  defined  as 


Px 

’  S*XX 

*  '^2*YY 

*  S*2Z 

^  S*YZ 

C  X 

6  XZ 

(25) 

Py 

’  "l^X 

^  ‘=2^Y 

^  ^3^22 

*  S^Y 

^  ^S^YZ 

c  v 

6^XZ 

(26) 

p. 

»  c  z 

♦  c.z,. 

♦  c  .z 

♦  c,z, , 

(27) 

2 

1  XX 

2  YY 

3  22 

4  XY 

5  YZ 

6  XZ 

The  determinant  of  the  Jacobian  transformation  matrix  Is  defined  as 


(2B) 


The  velocity  components  u,  v,  w  are  defined  as 

u  =  *  ^4'>.  * 

V  -  +  hj((.y  ♦  bg4>^)/0 

w  a  (h^^^  ♦  h^((>y  ♦  bg((>^;/D 


(29) 

(30) 

(31) 


A  second-order  local  coordinate  transformation  which  transforms  a  27-po1nt  cell 
from  the  physical  space  to  the  computational  space  (Fig.  1)  can  be  applied  to 
formulate  a  second-order  finite-difference  approximation  to  Eq.  (3)  as  des¬ 
cribed  ‘n  Ref.  3. 


(a)  Physical  space 
z 


(b)  Compuiaiional  space 


Figure  1.  Transformation  of  a  second-order  element. 


Equation  (3)  can  be  reduced  to  a  two-dimensional  equation: 

'l»XX  '  *  ‘4‘XY  •  '7*X  *  '8*Y  •  “■ 

With  c^,  c^.  c^,  c^,  and  Cg  redefined  as 

2  2  2  2  2 

-  [a^Xy  >  y‘)  -  U'^l/O^ 

2  2  2  2  2 

Cj  •  [a  (xj  ♦  Yx)  -  V  ]/0 


c^  *  -2[a^(Xj^XY  +  ^x^Y^  '  UV]/D^ 


C7  *  (X  Py  -  YyPx)/!^ 

's  -  ‘’x'x  - 

and  U,  V,  u,  V,  0,  Px,  and  py  redefined  as 


(32) 

(33) 

(34) 

(35) 

(36) 

(37) 


uyy  - 

"y 

(38) 

^^x  - 

"x 

(39) 

-  »*‘y''“ 

(40) 

^*X^Y 

U1 ) 

'x  =  ^’'XX  *  ^2Sr  ^  ^3*XY 


^l^X  ^  ^Z^YY  ^  '=3^XY 


('2) 

(O) 

(44) 


Equation  (32)  Is  consistent  with  the  two-dimensional  equation  derived  In  Ref, 
17.  The  detail  of  the  finite-difference  approximation  to  Eqs.  (3)  and  (32) 
can  be  found  In  Refs.  3  and  4.  * 


3 . 2  Finite-Volume  Hethod 

The  finite-volume  method  of  Jameson  and  Caughey  [5]  applies  the  mass  flux  con¬ 
servation  principle  to  formulate  the  finite-difference  approximation  to  the 
full-potential  equation.  As  shown  In  Fig.  2  for  a  two-dimensional  case,  four 
primary  cells,  1298,  2349.  9456  and  8967,  surrounding  the  control  point  9  are 
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Figure  2.  Primary  and  secondary  cells  In  the  f Inlte-vol jme  scheme. 


first  used  to  evaluate  the  flow  properties  ‘ncluding  the  velocity  components. 
Then  the  conservation  of  mass  flux  1s  applied  along  the  four  cell  boundaries 
of  the  secondary  cell  ABCO  wh'ch  1s  formed  by  the  centers  of  the  primary  cells. 
The  evaluation  of  the  mass  flux  on  these  four  secondary  cell  surfaces  Intro¬ 
duces  a  "lumping"  error  in  the  finite-difference  approximation  of  the  'ull- 
potentlal  equation  at  the  control  point  9,  In  addition,  the  finite-difference 
cDproximatlon  thus  constructec  results  In  uncoupling  of  even-  and  odd-rc'nt 
solutions.  The  remedy  to  this  problem  Is  to  add  proper  coupling  terms,  which 
essentially  shifts  the  evaluation  ot  mass  fluxes  from  the  centers  of  four  pri¬ 
mary  cells  to  the  control  point  9.  As  derived  In  tqs.  (;3)  ard  (34),  ‘re 
coefficients  c.|  and  c^  are  the  leading  terms  of  second  derivatives  of  the 
po'entlal  function  In  the  X  and  Y  direc-  tions,  respectively.  The  foUow'ng 
coupling  term  w'th  a  cross  derivative 


1  2  2  2  2 
»  2  (Xy  +  Yy)  <r  U 


compensates  for  the  lumping  error  In  computing  the  flux  balance  In  the  X-dlrec- 
tlon.  Similarly  the  following  coupling  term  with  a  cross  derivative 


12  2  2  2 
A,  .  j  la  (a,  •  y,)  -  V  )»„ 


compensates  for  the  lumping  error  In  computing  the  flux  balance  In  the 
Y-dIrectlon.  Similar  coupling  terms  can  be  constructed  for  three-dimensional 


cases . 
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with  these  coupling  terms,  the  finite-difference  approximation  of  the  finite- 
volume  method  Is  very  similar  to  that  of  the  finite-difference  method  des¬ 
cribed  In  Section  3.1  The  approximation  of  the  coupling  terms  In  the  finite- 
volume  method  Is  Identical  to  the  exact  coupling  term  In  the  finite-difference 
method  only  If  the  mesh  lines  are  perpendicular  to  each  other.  However,  since 
the  leading  terms  of  the  second  derivatives  of  the  two  method^,  are  Identical, 
the  relaxation  scheme  can  be  commonly  constructed  and  shared  by  both  methods 
while  allowing  different  ways  of  computing  the  residuals.  In  the  method  of 
Jameson  and  Caughey  [5],  the  coupling  terms  are  constructed  In  divergence  form. 
Therefore,  their  scheme  Is  fully  conservative.  Although  the  finite-difference 
method  has  an  exact  coupling  of  even-  and  odd-point  solutions,  the  method  Is 
rot  fully  conservative. 

Both  the  finite-difference  and  the  finite-volume  methods  are  Implemented  In 
the  present  study.  A  comparison  of  the  solutions  obtained  by  the  two  methods 
Is  given  In  Section  6.1. 

3.3  Artificial  viscosities  and  Shock-Point  Operators 

The  second  derivative  of  the  potential  function  In  the  st  earwise  direction, 

S,  Is  given  as 
2 

♦SS  *  J2  ^"^^x  *  ''Sy  *  *  ^^''♦xy  *  ^vw^^,  -  2uwo^^)  (47) 

Equation  (47)  can  be  rewritten  as 

^  ^  2VW^y2  ♦  (48) 

Q 

where  U,  V,  W  are  given  In  Eqs.  (13)  through  (15). 

The  directional  bias  of  supersonic  flows  carf  be  properly  simulated  by  perform¬ 
ing  an  upwind  differencing  or  adding  artificial  viscosities  In  the  approximate 
itreamwlse  direction.  If  Y  =  constant  lines  are  In  the  approximate  s  direc¬ 
tion,  the  principal  part  of  can  be  approximated  by 

^SS  ■  ^  '**XX 

q 

A  first-order  artificial  viscosity  can  be  expressed  by 
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XX. 


XX 


'1-1 


(50) 


where 


2 

y  .  max  (1  -  ^  ,  0)  (51) 

q 

»  • 

H  1s  then  added  to  the  finite-difference  representation  of  £q,  (3)  at  super¬ 
sonic  points.  At  shock  points,  I.e.,  the  first  downstream  subsonic  points 
after  the  shocks,  the  following  first-order  artificial  viscosity  Is  added 
with  p  controlling  the  nonconservative  differencing: 


H,  =  (P„  -  1)( 


uU^i> 


m 


Mx 
D  M- 


(52) 


m 

'w 


I 


If  p^  »  0,  the  quantity  yU  Is  conserved  along  Y  >  constant  lines,  1mp%/1ng 
that  the  added  artificial  viscosities  are  conserved  along  approximate  stream¬ 
lines.  If  p^  >  0,  a  numerical  mass  flux  Is  Introduced  at  shocks,  mod1fy1r*.g 

t'.e  locations  and  strengths  of  the  shocks.  The  effect  of  p  on  the  captu^^d 

m 

shocks  will  be  discussed  later.  A  second-order  artificial  viscosity  and  ‘nock- 
point  operator  can  be  expressed  as 


H  *  'AX)^  [ 


yU^(> 


XX. 


uU^<> 


XX 


'XX 


)l-2-  2  <■ 


yU^<> 


XX 


yU^^ 


)l-l  ^  (■ 


XX 


), 


(53) 


^  =  (Pm  - 


(AX)uU  ♦ 


XX. 


[(■ 


XX. 


yU^^ 


D  M-1 


-  (• 


XX 


)l-2^(Pm- 


'54) 


The  solution  Is  second-order  accurate  at  both  subsonic  and  supersonic  polints, 
and  first-order  accurate  at  shock  points.  Although  y  Is  a  rampf unction,  toth 
H  and  reduce  to  zero  as  the  mesh  size  goes  to  zero.  In  the  so-called 
quasi-conservative  schemes,  the  finite-difference  formulation,  described 
Section  3.1,  Is  applied,  and  only  the  differencing  of  artificial  vlscosltlees 
Is  In  divergence  form;  the  differencing  of  the  governing  potential  equatlm  is 
not.  A  second-order  fully  conservative  scheme  also  can  be  constructed  by 
Incorporating  H  and  Into  the  finite-volume  formulation  described  In  Set- 
tlon  3.2. 


similarly,  a  nonconservative  scheme  can  be  developed  by  adding  the  following 
artificial  viscosity  at  supetsonlc  points 


H  •  (  0 

to  make  the  differencing  upwind  as  In  the  original  Murman  and  Cole  scheme  [18], 
while  no  artificial  viscosity  Is  added  at  the  shock  point. 
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4.0  CLCBSCH  TRANSFORMATION 
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The  Clebsch  transformation  has  been  extensively  employed  In  the  calculation  of 
three-dimensional  rotational  flows  in  turbomachines  up  to  the  subcritical  flow 
range  using  mostly  analytical  approaches.  The  method  has  ' ?en  successfully 
developed  for  computing  shear  flows  [19,20],  wake  flows  [21,22,23],  and  non- 
axl symmetric  Inlet  flows  [24], 

In  the  Clebsch  formulation  of  steady  rotational  flows,  the  velocity  vector  v  Is 
decomposed  Into  a  potential  part  and  rotational  parts,  written  In  terms  of 
scalar  functions  [19]: 


V  .  (56) 

Hence,  by  definition,  the  vo'^tlcUy  ector  Q  Is  nonzero: 

Q  =  VxV  .  I  Vff  X  VX  (57) 

~  "  n  n  n 

To  determine  the  flowfleld,  the  Clebsch  variables  ♦,  a  and  X  must  be  chosen 

n  n 

so  that  the  equations  of  motion  are  satisfied.  In  general,  each  pair  of  a  and 

n 

X^  can  be  considered  to  represent  the  vorticity  field  generated  by  various 
shocks  or  the  trailing  vorticity  downstream  of  lifting  bodies  such  as  wings  or 
propeller  blades. 


For  steady  flow. 


V  •  (pl)  -  0 


while  the  momentum  equation,  written  In  Lamb's  form  for  Isoenergetic  flow  In 
tr.e  absence  of  body  forces.  Is 


2 

V  X  Q  »  -  —  7s  (59) 

-  -  T  , 

where  p  Is  the  density,  s  the  entropy,  a  the  speed  of  sound,  and  t  the 
ratio  of  the  specific  heats  of  the  flow. 


In  the  present  approach,  the  potential  part  In  Eq.  (56)  Is  determined  by  solv¬ 
ing  the  full-potential  equation  given  1n  the  previous  section,  while  the  rota¬ 
tional  parts  In  Eq.  (56)  are  chosen  to  satisfy  the  momentum  equation,  Eq.  (59) 
The  entire  flowfleld  can  then  be  computed  by  solving  the  governing  coupled 
equations  of  these  Clebsch  variables  iteratively. 
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When  shocus  appear  In  the  flowfleld,  entropy  Increases  across  the  shocks.  The 
rotational  parts  of  the  velocity  vector  defined  In  Eq.  (56)  are  determined 
from  the  momentum  equation.  A  general  solution  to  Eq.  (59)  Is  given  by: 

Q  »  Q.  ♦  n  (60) 

~H  -p 

Here  Is  the  homogeneous  solution  of  Eq.  (59),  namely: 

V  X  -  0  (61) 

and  Q  Is  the  particular  solution  of  Eq.  (59),  namely: 

~P 

V  X  Qp  .  -  ^  7s  (62) 

The  particular  solution  fl  represents  the  vorticity  component  which  Is  not 
parallel  to  the  velocity  vecto'. 

Following  arguments  similar  to  those  1n  Refs.  19  to  24,  the  vorticity  compon¬ 
ent  associated  with  the  entropy  jump  across  shocks  can  be  expressed  as: 

Q  .  7t  X  7s  (63; 

~P 

where  s  Is  the  entropy  field  and  t  Is  a  Clebsch  variable.  Substituting  Eq. 

(62)  Into  Eq.  (61),  and  considering  that  entropy  Is  convected  along  stream¬ 
lines  behind  shoes, 

V  •  7s  =  0  (54) 

It  can  be  shown  that  the  goveining  equation  for  t  Is 

2 

V  •  7t  =  —  (FF) 

"  T 

The  Clebsch  variable  t  Is  similar  to  the  Oarwin-Llghthlll-Hawthorne  drift 
function  [19].  The  variation  of  x  from  streamline  to  streamline  Is  directly 
connected  to  the  stretching  and  tipping  of  the  vortex  filaments  associated 
with  the  entropy  variation  In  the  flowfleld.^ 

The  homogeneous  solution  represents  the  vorticity  component  which  Is  parallel 
to  the  velocity  vector.  For  the  problem  considered  here,  Q  Is  the  trailing 
vorticity  shed  behind  a  wing.  In  the  airfoil  case.  Is  Identically  zero. 

It  can  be  shown  that  the  homogeneous  solution  Is 

,  4(E)7E  X  7r  (66) 
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where  5  5  y  -  f(x,z)  «  0  defines  the  wake  surface  or  the  location  of  the 
trailing  vortex  sheet,  4(5)  Is  the  Oirac  delta  function,  and  f  Is  the 
circulation  of  the  vortex  filaments  shed  by  the  wing.  Note  that  according  to 
Eg.  (66),  the  trailing  vorticity  lies  In  the  wake  surface  and  Is  zero  every¬ 
where  except  on  the  wake  surface.  Substituting  Eq.  (66)  Into  Eq.  (61),  along 
with  the  use  of  the  wake  boundary  condition,  namely: 

‘  H  •  0  (67) 

the  condition  for  zero  pressure  Jump  across  the  wake.  Including  the  trailing 
edge.  Is  obtained 

V  •  7r  -  0  (68) 

Here  v  Is  the  mean  velocity  between  the  wake  upper  and  lower  sides. 


Given  the  velocity  field  from  the  previous  Iteration  step,  the  Clebsch  vari¬ 
ables  s,  T.  I  and  r  can  be  updated  by  solving  Eqs.  (64).  (65),  (67)  and 
(68),  respectively. 


As  mentioned  earlier,  the  full-potential  equation  Is  employed  to  determine  the 
potent'jl  part  of  the  velocity  vector.  Using  the  results  derived  In  Section 
2.1,  the  velocity  vector  can  be  written  as: 

V  »  -  sVt  ♦  H(5;7r  (69) 

where  H(E)  Is  the  step  function.  Substituting  Eq.  (69)  Into  Eq.  (58)  /lelds 


V  •  {p[7<>  -  s7t  ♦  H(E)7r]}  .  0 


(70) 


Given  the  rotational  parts  from  the  previous  Iteration  step,  the  potential 
part  can  be  updated  by  solving  Eq.  (70)  in  t^e  following  form: 


V  .  -  7  .  (p{-s7t  ♦  H(5)Vr)}  (71) 

The  density  p  can  be  related  to  the  local  flow  properties  by  means  of  the 
isentropic  relation  and  the  energy  equation: 


and 


(■>2) 
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2  1  T  -  1  ,1  2. 

a  ■  -T  ♦  -"-J—  (1  -  q  ) 

m 

where  q  Is  the  magnitude  of  the  velocity  vector,  and  Is  the  freestream 
Mach  number. 
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5.0  NUMERICAL  PROCEDURE  AND  RESULTS  FOR  TWO-DIMENSIONAL  FLOWS 


The  governing  equations  and  boundary  condition,  the  numerical  procedures  to 
solve  them  and  the  results  for  two-dimensional  flows  are  presented  In  this 
section. 

»  ‘ 

5 . 1  Governing  Equations 

The  finite-difference  and  finite-volume  methods  can  now  be  modified  to  Include 
the  effects  of  the  total  pressure  loss  across  shocks  and  the  vorticity  down¬ 
stream  of  shocks.  For  flow  around  an  airfoil,  there  Is  no  trailing  wake  down¬ 
stream  of  the  trailing  edge.  From  Eqs.  (40),  (41)  and  (69),  and  applying  the 
same  coordinate  transformation  as  1n  Section  2.0,  the  velocity  components  can 
be  expressed  as 


V  -  [(*x'^Y  '  *y'^X^  *  ^^*x''y  ■ 


(74) 


and  the  governing  equations  can  be  expressed  as 

^I'^XX  *■  *'2^YY  *  ''4'^XY  *  ’'T^X  *  ^8^Y  "  ’•(P^f)/p  (75) 

with  the  following  two  equations  for  solving  the  Clebsch  variables  s  and  r: 

USX  ♦  VSy  .  0  (76) 

2 

♦  Vr^  -  D  ^  (77) 

A  multigrid  line-relaxation  scheme,  originally  developed  for  transonic  full- 
potential  methods  [24-28]  Is  applied  to  solve  £q.  (75)  for  the  potential 
function  ♦,  along  with  the  new  expression,  Eq.  (74),  to  coeipute  the  velocity 
components.  The  Clebsch  variables  can  be  solved  analytically  under  certain 
approximations,  or  numerically  by  using  the  Lax-Wendroff  scheme  [29]. 

5.2  Approximate  Euler  Clebsch  Approach 

The  flowfleld  around  an  airfoil  can  be  divided  Into  Irrotatlonal  and  rotational 
regions  as  shown  In  Fig.  3.  The  flowfleld  downstream  of  the  shock  Is  rota¬ 
tional.  Across  a  shock,  there  is  an  entropy  Jump  which  c»n  be  estimated  from 
the  Rankine-Hugonlot  shock  relation  1f  the  Mach  number  upstream  of  the  shock 
Is  known.  The  entropy  s  Is  then  convected  downstream  along  the  streamlines  as 
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Figure  3.  Irrotatlonal  and  rotational  Fiowfleld  about  an  airfoil. 

described  In  £q.  (76).  In  the  present  method,  a  C-.nesh  which  conforms  to  :.'.e 
airfoil  surface  Is  applied.  Therefore,  It  Is  a  fairly  good  approximation  that 
the  entropy  Is  constant  along  the  mesh  lines  Y  =  constant  downstream  of  the 
shocx . 


A  second  approximation  can  be  made  In  order  to  solve  for  the  Clebsch  variable 
T  analytically.  As  described  In  Section  4.0,  t  represents  a  material 
coordinate  s.-face  which  stretches  according  to  the  local  flow  velocity,  as 
described  In  Eg.  (77).  If  the  velocity  is  assumed  to  be  freestream  ve'jcity, 
and  the  speed  of  sound  Is  also  assumed  to  be  freestream  value.  Eg.  (77)  can  be 
'ewritten  as 


,  3t  .  5t  1 

(ccsa)  -  *  (sina)  -  =  -j 


(78) 


where  a  Is  the  angle  of  attack. 


In  order  to  solve  Eg.  (78),  boundary  conditions  need  to  be  prescribed  for  t 
at  the  shock  front.  Since  the  potential  function  ^  Is  taken  to  be  continu¬ 
ous  across  the  shock,  the  conservation  of  tangential  momentum  across  the  shocx 
front  requires  t  to  be  a  constant  everywhere  on  the  shock  surface.  However, 
since  the  Clebsch  variable  i  Itself  Is  not  used  In  the  present  formulation 
to  compute  the  velocity,  but  rather  Its  gradient,  the  boundary  condition  Is 
enforced  by  requiring  that  the  Jump  In  the  velocity  vector  across  the  shock  Is 
normal  to  the  shock  surface.  Since  all  the  Jump  In  the  velocity  across  the 
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shock  Is  contained  In  the  rotational  part,  this  condition  1s  satisfied  by  the 
fol lowing: 

7t  X  n  «  0  (7« 

where  n  =  (^x^®x  *  local  shock  surface  unit  normal.  Solv¬ 

ing  tqs.  (78)  and  (79)  simultaneously  for  the  gradient  of  i  gives.  In  the 
two-dimensional  case. 


COSO  ♦  n^  sino) 

Hence,  In  the  2-0  Euler  correction  method,  the  velocity  vector  reduces  to 


V  .  -  — -  n  (( 

COSO  ♦  n^  sino) 

•  X  y 

since  Vt  Is  constant  along  Y  =  constant  lines,  the  source  term  on  the  rlght- 

har,  side  of  Eq.  (75)  Is  small  everywhere  except  at  shock  points  where  dis¬ 
continuities  In  s  and  p  exist. 

To  obtain  the  flowfleld,  the  following  Iterative  procedure  Is  employed; 

1.  Set  the  Initial  solution  for  ♦  to  be  the  freestream  condition,  and  set 
s  and  T  to  be  zero. 

2  Solve  Eq.-  (75)  for  using  the  multigrid  relaxation  met'oc. 

3.  Obtain  the  shock  surface  normal  vector  and  the  entropy  Jump  across  the 
shocks  from  the  Rankine-Hugonlot  relation  when  supersonic  pockets  (or 
shocks)  start  to  appear  In  the  flowfleld. 

4.  Update  the  rotational  part  In  Eq.  (81). 

5.  Repeat  Steps  2  to  5  until  solutions  converge. 


5.3  Exact  Euler-Clebsch  Approach 

The  governing  equations  for  s  and  t,  Eqs.  (76)  and  (77),  are  of  hyperbolic 
type.  They  can  be  solved  using  the  second-order  Lax-Wendroff  explicit  scheme. 
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Consider  the  following  general  equation; 


where  F 


=  (s.t)  and  G  r  (0,  a  /y) 


(82) 


The  solution  procedure  to  solve  the  above  equation  consists  of  prescribing  the 
Initial  condition  for  F  at  the  shock  front,  and  solving  Eq.  (82)  by  marching 
downstream  In  the  X  direction.  During  space  marching,  the  solution  of  F  at 
the  1  station  Is  given,  and  the  solution  at  the  (1  ♦  1)  station  Is  computed  In 
two  steps  [29). 


Step  1 


1  ''l 


.1 


.1, 


1  S 


pUl/2  i  ,f1  ^1,  ^  ^  ... 

^J*l/2  *  2  ^^J*l  *  ■  2  *  2  U 


fUV:  i  ,p1  p1  ,  1  !i  ,.1  .1  ,  1  ^ 

^J-1/2  *  2  ^J-1^  “  2  ‘  ^J-1^  *  2 


(93) 


(84) 


Step  2 


V  V  G  G 

.  1  (!i  ,  .  p1*l/2  ,  (!i  , 

J  2  'u^  j+1/2  j-1/2^  'u^ 


35) 


The  subscripts  1  and  2  denote  the  cell-centered  va'ues  of  the  cells  defined  in 
Fig.  4.  In  this  scheme,  the  Von  Neumann  stability  requ^'-es 


I-' 


(86) 


In  the  flowfleld  where  the  scheme  applies,  the  surface  normal  velocity  compon¬ 
ent  Is  usually  small  compared  with  the  streamwise  component. 

The  coupling  between  the  potential  and  rotational  solution  procedures  Is  sim¬ 
ilar  to  the  one  described  In  Section  5.1,  except  In  Step  4  where  the  numerical 
solutions  are  computed  for  s  and  t  Instead  of  using  the  analytical  solutions. 

5.4  Kutta  Condition 

The  Kutta  condition  requires  that  the  static  pressures  at  the  upper  and  lower 
trailing  edges  be  matched.  Furthermore,  It  also  requires  that  there  Is  no 
pressure  Jump  across  the  streamline  emanating  from  the  trailing  edge.  In  the 
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Figure  4.  Lax-Wendrof f  scheme. 


full-potential  method,  flow  Is  assumed  lo  be  Isentroplc  everywhere.  Therefore 
the  Kutta  condition  can  be  satisfied  by  assuming  the  potential  Jump  across  the 
trailing  edge  Is  constant  along  the  streamline,  or  the  mesh  line,  emanating 
from  the  trailing  edge.  This  ensures  that  the  magnitudes  of  the  velocities  on 
both  sides  of  the  streamline  are  equal,  as  are  tne  pressures  on  both  sides  of 
the  streamline.  In  the  present  method  If  shocks  appear  on  either  upper  or 
lower  or  both  surfaces,  the  total  pressures  of  the  streamlines  Just  above  and 
below  the  trailing  edge  are  different.  Hence,  in  order  to  have  equal  pressure 
across  the  dl'viding  streamline,  the  velocity  must  be  discontinuous  across  It. 

In  this  section  the  velocity  discontinuity  across  this  dividing  streamline  will 
be  shown  to  come  from  the  rotational  term  sVt  so  that  the  condition  for  ^ 
across  this  dividing  streamline  Is  the  same  as  In  the  full-potential  method. 

To  simplify  the  derivation,  consider  the  case  where  there  Is  only  one  shock  on 
the  airfoil  upper  surface.  The  static  pressure  can  be  written  In  terms  of 
local  velocity  and”  stagnation  pressure: 


^tl  Y  -  1  2  2 


(87) 


Y/Y-1 


(88) 
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Equating  Eqs.  (87)  and  (88) 


['  •  -  =i?)i 


[1  ♦  0  -  q2)  .  (q  .  i)M2q^&q 


V  1  ?  5  - 

-  M;Aq‘^]  (89) 

2 

where  fiq  =  q^  -  q^ .  Neglecting  the  term  Aq  In  Eq.  (89),  an  explicit  expres¬ 
sion  for  the  velocity  Jump  across  the  dividing  streamline  can  be  found. 


1  1  ^t?  Y  -  1  2  2 

Aq  =  y  -y  [1  -  (—)][!  ♦  (1  -  q()]  (90) 

‘’l  ^tl  ' 

09 

In  practice,  the  stagnation  pressure  loss  across  the  shock  Is  gcierally  small 

even  for  a  strong  shock.  For  example,  a  shock  Mach  number  of  1.7  results  In 

15%  loss  In  stagnation  pressure.  Since  Aq  Is  proportional  to  the  stagnation 

2 

pressure  loss,  the  assumption  that  the  Aq  term  can  be  neglected  In  deriv¬ 
ing  Eq.  (90)  Is  valid. 


From  Eq.  (69),  the  velocity  Jump  across  the  dividing  streamline  is  given  by 

Aq  »  (s  .5.) 

where  2  Is  the  distance  measured  along  the  dividing  streamline. 


For  1 soenergeilc  flow,  the  entrooy  1s  related  to  the  stagnation  pressure  by 


s 


1  -In 


't2, 


'tl 


=  1 


XI 

p 


(92) 


where,  as  before,  second-order  terms  In  the  stagnation  pressure  loss  are  neg¬ 
lected  . 


From  the  governing  equation  of  the  Clebsch  variable  i,  Eq.  (65), 


3t 

dl 


-  —  [1 
<^1  ym2  ^ 


(1  -  q^] 


(93) 
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The  numerical  solutions  obtained  by  the  present  method  for  airfoils  are  pre¬ 
sented  In  this  section.  Figure  6  shows  the  comparison  between  the  solutions 
obtained  for  the  NACA-0012  airfoil  at  M  »  0.800  and  a  =  1.25*  from  the  two 
Euler-Clebsch  approaches  and  the  time-marching  Euler  method  [!,]•  The  fully 
conservative  full-potential  solution  predicts  a  shock  almost  at  the  trailing 
edge.  The  source  term  In  Eq.  (75)  in  the  present  method  accounts  for  correc 
tions  to  the  total  density  and  the  vorticity  which  are  neglected  In  the  full 
potential  solutions.  Solutions  obtained  by  using  the  approximate  and  the 
exact  Euler-Clebsch  approach  predict  weaker  and  more  upstream  shocks  which 
agree  better  with  the  time-marching  Euler  solution.  A  study  of  the  order  of 
magnitude  on  Eqs.  (74),  (75)  and  (76)  reveals  that  the  first-order  correct'o 
comes  from  the  source  term  of  Eq.  (75)  at  shock  points  where  entropy  Jump 
occurs,  and  that  the  second-order  correction  comes  from  the  vorticity  effect 
which  modifies  the  velocity  distribution,  Eqs.  (74)  and  (75),  and  source 


-1.0 


o  0.0 

u 


EULER  (AGARO-AR-211 ) 

APPROXIMATE  EULER  CLEBSCH 

EXACT  EULER  CLEBSCH 

FULLY  CONSERVATIVE  FULL  POTENTIAL 


Figure  6.  Comparison  of  solutions  for  NACA-0012  airfoil  at  *  0.8  anq 
a  .  1.25*. 
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term  of  Eg.  (76)  at  points  downstream  of  shock  points.  In  the  exact  Euler- 
Clebsch  approach,  the  rotational  part  Is  evaluated  more  accurately,  and  the 
solution  seems  to  agree  better  with  the  time-marching  Euler  solution  than  the 
approximate  Euler-Clebsch  solution.  Note  that  the  shocks  captured  by  the 
Euler-Clebsch  solutions  are  not  as  sharp  as  those  obtained  by  the  time-marching 
Euler  solution,  especially  on  the  lower  surface.  This  Is  probably  because  the 
time-marching  Euler  computation  employs  more  points  on  the  airfoil  surfaces 
than  the  Euler-Clebsch  solutions  (192  points  vs  162  points). 


Figure  7  presents  the  solutions  obtained  for  the  RAE-2822  airfoil  at  M  . 

dD 

0.75  and  a  «  3*.  This  test  case  was  chosen  In  order  to  evaluate  the  accu¬ 
racy  of  the  Euler-Clebsch  method  In  the  presence  of  a  strong  shock.  Figure  3 
Illustrates  the  stagnation  pressure  contour  obtained  using  an  In-house  version 


EULER-CLEBSCH  METHOD 


TIME-MARCHING  EULER  METHOD 


Figure  8.  Comparison  of  total  pressure  contours. 

of  the  Jameson  DFLO-53  t1me-march1ng  Euler  method  [30,31]  and  the  exact  Euler- 
Clebsch  solution.  This  figure  clearly  shows  that  the  convectUe  behavior  of 
the  stagnation  pressure  along  streamlines  Is  preserved  In  the  present  Euler- 
Clebsch  method.  The  time-marching  Euler  method,  however,  has  difficulty 
modeling  this  Inviscid  characteristic  because  of  the  presence  of  numerical 
dissipation  terms. 

The  present  exact  Euler-Clebsch  method,  however,  has  difficulties  In  giving  a 
converged  solution  when  the  shock  becomes  too  strong  or  the  angle  of  attack 
becomes  too  large  due  to  the  stability  criterion  of  Eg.  (86).  More  stable 
schemes  such  as  the  Crank-Nichol sen  Implicit  scheme  are  currently  under 
Investigation. 
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6.0  NUMERICAL  PROCEDURE  AND  RESULTS  FOR  THREE-DIMENSIONAL  FLOWS 

The  numerical  procedure  to  solve  the  governing  equation  with  boundary  condi¬ 
tions  for  the  three-dimensional  flows  Is  discussed  In  Section  6.1.  The  numer¬ 
ical  results  obtained  for  wings  and  wing/bodies  are  presented  In  the  subsequent 
sections.  A  comparison  of  finite-difference  and  finite-volume  solutions  ls 
presented  In  Section  6.2.  Both  full-potential  and  Euler-Clebsch  solutions 
obtained  for  wings  are  presented  In  Section  6.3,  and  the  solutions  obtained 
for  wing/bodies  are  presented  In  Section  6.4. 

6 . 1  Governing  Equations 

In  the  three-d 1me''s Iona  1  case,  the  approach  Is  similar  to  the  two-dimensional 
case,  except  for  some  additional  assumptions  and  modifications.  Only  the 
approximate  Euler-Clebsch  approach  Is  applied  here.  First,  the  waxe  or  the 
trailing  vortex  sheet  defined  by  E  In  Eq.  (67)  Is  taken  to  be  the  grid  surf¬ 
ace  emanating  from  the  wing  trailing  edge.  Furthermore,  the  trailing  vorticity 
Is  assumed  to  be  convected  along  the  grid  lines  leaving  the  trailing  edge. 

This  Is  done  by  setting  the  Jump  In  ^  across  the  wake  at  each  spanwise  sta¬ 
tion  to  be  constant  along  the  mesh  lines  downstream  of  the  trailing  edge. 

Hence,  the  trailing  vorticity  Is  the  homogeneous  solution  defined  In  Eq.  (66). 
This  treatment  of  the  wake  Is  the  same  as  that  applied  In  most  potent'al 
aporoaches . 


Second,  for  swept  wings,  oblique  shocks  are  usually  found  on  most  parts  of  the 
wing  surfaces.  In  some  cases,  1t  1s  possible  to  encounter  oblique  shocks  where 
the  total  Mach  numbers  both  upstream  and  downstream  of  these  shocks  are  super¬ 
sonic.  Although  this  type  of  oblique  shock  does  appear  In  many  transonic 
applications.  Its  shock  strength  Is  generally  weak;  therefore,  no  entropy  cor- 

ft 

rectlon  Is  given  to  this  type  of  shock.  In  the  absence  of  yaw,  Eqs.  (79)-(81) 

still  hold  In  3-0  where  the  local  shock  surface  unit  normal  vector  Is  defined 

asn*(n)e  ^Inle  ♦■(n)e. 

'  x'  X  y  y  '  z  z 

The  modified  full -potential  equation  can  now  be  derived  from  Eqs.  (3)  and  (50) 
as 


KXX 


c^4»yy  ♦  ^ 

2 

♦  Cgi>^  »  Da  V  •  (psVt)/p 
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^7^X 


(95) 


wUh  the  velocity  components,  defined  as 


u  •  ♦  ^4'>Y  "  ^7^2^  *  *  ^4‘^Y  * 

V  •  "■  ^8*^2^  "■  ^  ^5^Y  *  ^8^2^^^° 

w  -  [(^3^  "  ^♦y  *  V2^  *  "  S"2^^f.° 

In  the  approximate  Euler-Clebsch  approach,  s  Is  assumed  to  be  constant  along 
Y  .  constant  lines.  Following  the  same  argument  as  In  the  two-dimensional 
case,  the  entropy  s  can  be  computed  by  the  Rankine-Hugonlot  relation,  and  the 
gradient  of  t  can  be  approximated  by  Eq.  (80). 

In  the  downstream  farfleld,  the  Neumann  boundary  condition  1s  employed.  The 
Implementation  of  this  type  of  boundary  condition  requires  the  knowledge  of 
the  far  downstream  velocity,  which  Is  not  known  a  priori  and  Is  a  function  of 
the  shock  strength.  In  the  present  study,  the  far  downstream  velocity,  denoted 
by  V43 .  Is  taken  to  be 

=  (cosa)e  ♦  (slna)e  -  — r - - -  n  (97) 

^  yW^  (f'  COSO  ♦  n  sintt) 

•  X  y 

Conceptually,  the  last  term  In  Eq.  (97)  Is  related  to  the  wave  drag,  and  rep¬ 
resents  the  velocity  deficit  due  to  the  loss  In  the  fluid  stagnation  pressure 
across  the  shock. 

6  2  Comparison  of  F*n1te-01fference  and  Finite-Volume  Solutions 
As  described  in  Section  3.2,  two  major  differences  In  the  finite-difference 
and  the  f 1 ni te- vol ume  methods  come  from  the  formulation  of  the  odd-  and  even- 
point  coupling  terms,  and  the  calculation  of  the  density  term.  The  coupling 
term  In  the  finite-difference  method  Is  derived  directly  from  the  coordinate 
transformation,  while  the  coupling  term  In  tfe  f 1 n1 te-vol ume  method  accounts 
for  the  leading  terms  only.  The  explicit  calculation  of  the  density  term  In 
the  f 1 nl te- VO  1 ume  method  makes  the  scheme  conservative,  but  requires  more 
computational  work  because  of  the  need  to  compute  exponential  quantities. 

Figures  9  and  10  present  the  solutions  obtained  by  the  finite-difference  and 
f 1 nl te- VO  1 ume  methods.  The  surface  boundary  condition  In  the  finite-difference 
method  Is  enforced  exactly  as  It  Is  In  the  finite-volume  method,  while  differ¬ 
ent  formulations  are  used  In  evaluating  the  residuals.  Both  solutions  are 
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figure  9.  Comparison  of  f Inite  difference  and  finite  volume  soluti 


Obtained  using  the  multigrid  scheme.  Figure  9  presents  two  solutions  obtained 
for  a  NASA  swept  wing  [1]  at  «  0.833  and  o  -  1.75*.  The  finite-difference 
and  finite-volume  solutions  agree  well  except  for  minor  differences  In  the 
after-shock  flow  reexpansion  near  the  wing  tip.  Figure  10  presents  a  compari¬ 
son  of  three  computed  solutions  obtained  by  the  present  method  and  FLO-22  [32] 
for  an  0NERA-M6  wing  at  M  »  0.84  and  a  =  3.06®.  The  nonconservative  scheme  Is 

OB 

applied  In  the  present  method  In  order  to  make  a  fair  comparison  with  the 
FLO-22  solution  which  Is  also  nonconservative.  Both  the  present  finite-volume 
and  finite-difference  solutions  are  obtained  using  160  x  24  x  32  mesh  points, 
while  the  FLO-22  solution  Is  obtained  using  192  x  24  x  32  points.  The  FLO-22 
method  applies  an  extrapolated  relaxation  scheme  to  accelerate  the  solution 
convergence,  while  the  present  method  applies  a  multigrid  scheme.  In  general, 
the  present  solutions  converge  better  than  t'e  FLO-22  solutions.  The  agreement 
between  the  three  solutions  Is  generally  good,  except  for  a  small  discrepancy 
In  the  prediction  of  shock  locations. 

Figure  11  shows  a  similar  comparison  of  the  solutions  obtained  for  a  Douglas 
wind-tunnel  model  wing,  LB-488.  The  LB-488  wing  has  a  significant  aft  loading. 
The  present  solutions  predict  an  almost  flat  pressure  plateau  near  the  wing 
root,  while  FLO-22  predicts  a  slight  oscillation.  Most  Interestingly,  both  the 
present  finite-volume  and  finite-difference  solutions  predict  a  significant 
preshock  reexpansion  between  the  65  to  901C  semispan  locations,  while  the  FlO-22 
solution  does  not.  This  preshock  reexpansion  was  first  observed  In  the  finite- 
volume  solution,  and  It  was  uncertain  then  whether  this  reexpansion  was  due  to 
numerical  oscillation  or  possibly  due  to  the  special  coupling  terms  uescribed 
1n  Section  3.2.  The  present  finite-difference  solution  obtained  with  more 
exact  coupling  terms  still  shows  the  preshock  reexpansion.  Further  Investiga¬ 
tion  with  the  use  of  finer  meshes  should  be  carried  out  In  order  to  understand 
the  discrepancy  between  the  present  solution!  and  the  FLO-22  solution. 

6 . 3  Euler-Clebsch  Solutions 

The  Euler-Clebsch  solutions  are  obtained  using  the  analytical  approach  des¬ 
cribed  In  Section  6.1.  Figures  12  to  14  present  the  solutions  obtained  for 
the  NASA  swept  wing.  A  comparison  of  the  pressure  distributions  computed  at 
several  spanwise  stations  using  Jameson's  time-marching  Euler  method,  the  pres¬ 
ent  Euler-Clebsch  method,  and  the  nonconservative  full-potential  method  Is 
shown  In  Fig.  12.  This  figure  clearly  demonstrates  that  the  Euler-Clebsch 
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solution  ls  In  good  agreement  with  the  time-marching  Euler  solution,  except 
near  the  wing-tip  region.  On  the  other  hand,  the  nonconservative  full- 
potential  method  tends  to  predict  earlier  shock  formation  and  weaker  shock 
strength  compared  to  the  time-marching  Euler  method.  Moreover,  the  after-shock 
reexpansion  phenomenon  appears  In  the  nonconservative  full-potential  solution 
when  the  shock  becomes  strong,  for  example,  at  the  40%  and  60%  span  locations. 

Near  the  wing-tip  region,  both  the  full-potential  and  the  Euler-Clebsch  methods 
predict  a  higher  suction  peak  than  the  time-marching  Euler  method.  Moreover, 
large  after-shock  reexpansions  are  found  In  the  solutions  of  the  two  former 
methods.  Figure  13  gives  comparisons  of  sectional  C.,  and  C  between  the 
results  obtained  using  these  three  methods.  Finally,  Figure  14  shows  the  con¬ 
vergence  histories  of  average  residuals  and  percent  of  supersonic  points  In  the 
computational  domain  as  a  function  of  work  unit  by  the  full-potential  and 
Euler-Clebsch  methods.  This  figure  clearly  demonstrates  that  the  Euler-Clebsch 
method  has  the  same  convergence  rate  as  the  full-potential  method  with  about 
20%  more  computational  time  per  work  unit.  Moreover,  the  computer  storage 
requirements  In  these  two  methods  are  approximately  the  same.  On  the  other 
hand,  the  t1me-march1ng  Euler  method  Is  estimated  to  require  nearly  one  order 
of  magnitude  more  computational  time  than  the  ful l-pctent1al  method. 

•t  zero  angle  of  attack,  the  0NERA-M6  wing  Is  nonllfting  and  hence  does  not 
nave  a  trailing  vertex  sheet.  Inaccuracies  associated  with  the  approximate 
modeling  of  tne  trailing  vortex  sheet  can  therefore  be  Isolated  In  this  test 
case.  Figure  15  Illustrates  a  comparison  of  the  pressure  distributions 
obtained  at  three  spanwise  locations  using  the  Euler-Clebsch  method,  the  non- 
conservative  full-potential  method,  and  the  Onera/Matra  time-marching  Euler 
method.  The  Onera/Matra  solution  presented  In  the  AGARO  report  [1]  was  chosen 
for  comparison  since  It  lies  between  the  oth%r  AGARO  solutions  and  has  tabu¬ 
lated  pressure  data.  In  general  the  Euler-Clebsch  method  predicts  shock  loca¬ 
tion  and  shock  strength  better  than  the  full-potential  method.  The  shock 
locations  predicted  by  the  Euler-Clebsch  method  are  consistently  slightly  down¬ 
stream  of  those  predicted  by  the  time-marching  Euler  solution,  except  near  the 
wing-tip  region.  For  example,  at  the  80%  span  station,  the  Euler-Clebsch  solu¬ 
tion  predicts  the  shock  location  slightly  upstream  of  that  predicted  by  the 
time-marching  Euler  solution.  In  addition,  the  study  also  Indicates  that  the 
Euler-Clebsch  method  predicts  a  suction  level  In  front  of  the  shock  slightly 
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Figure  14.  Comparison  of  convergence  history  for  the  solutions  obtained  for 
NASA  swept  wing  at  =  0.833  and  a  *  1.75®. 


higher  than  that  predicted  by  the  time-marching  Euler  solution,  hence  the  shock 
Is  stronger. 


All  of  these  discrepancies  are  very  similar  to  those  found  In  the  NASA  swept 
wing  case,  even  though  the  trailing  vortex  sheet  is  absent  1n  this  example. 
Hence,  It  1s  believed  that  the  disagreements  around  the  wing-tip  region  In  the 
results  of  these  methods  may  be  associated  with  the  different  ways  In  which  the 
wing-tip  boundary  conditions  are  handled.  However,  due  to  the  limited  data 
available  In  the  literature  and  the  premature  status  of  the  time-marching  Euler 
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code  at  the  present  time,  no  definite  conclusions  can  be  drawn  from  these 
results  as  to  the  accuracy  of  the  present  method  near  the  wing  tip  region. 

6.4  W1nq/6odv  Solutions 

An  example  of  a  wing/body  flowfleld  calculation  Is  demonstrated  for  an  F-14 
model  [33].  The  Input  geometry  Is  shown  In  Fig.  16.  Near  thewing  root,  there 
Is  a  leading-edge  break,  where  the  leading-edge  sweep  angle  changed  from  68*  to 
22“.  The  Inlet  entrance  Is  modeled  by  a  flat  surface  where  the  no-flux  condi¬ 
tion  Is  applied.  The  grid  generation  method  of  Chen,  Vassberg  and  Peavey  [28] 
Is  applied  to  generate  a  C-H-H  grid,  as  shown  In  Figs.  17  and  18.  Figure  17 
shows  a  grid  distribution  on  the  wing  and  fuselage  surfaces,  while  Fig.  18 
shows  a  typical  fuselage  cross-sectional  grid  distribution  at  nearly  midwing. 
The  solutions  obtained  for  wing-alone  and  wing/body  cases  are  presented  In 
Figs.  19  and  20  for  M  =  0.800  and  a  =  4*,  and  M  »  C.900  and  a  »  2*,  respec- 

CO  CO 

tively,  using  the  nonconservative  scheme.  The  results  show  that  the  fuselage 
Introduces  a  significant  upwash  effect  on  the  pressure  distribution  on  the  wing 
surface.  This  effect  Is  more  pronounced  near  the  wing  '•oot  than  near  the  wing 
tip.  The  slight  oscillation  In  the  wing/body  solutions  near  the  leading  edge 


i., 


Figure  16.  An  Input  geometry  for  F-14  wing/body. 
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figure  17.  Surface  grid  distribution  on  the  F-U  wing/body 
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Is  probably  due  to  the  mesh  distortion  caused  by  an  unsmooth  variation  of  the 
fuselage  cross-sectional  area. 

Comparisons  of  solutions  obtained  for  the  r-14  wing/body  using  the  fully  con¬ 
servative  ful 1 -potential  method  and  the  present  Euler-Clebsch  method  are  pre¬ 
sented  In  F^gs.  21  and  22  for  M  .  0.85  and  a  .  4*,  and  «  •  0.9  and  a  =  2*. 

m  m 

respectively.  Figure  21  shows  that  the  f u1 ly-conservatl ve  full-potential  solu¬ 
tion  predicts  stronger  and  more  downstream  sfiocks  than  the  Euler-Clebsch  solu¬ 
tion.  An  exceptionally  high  suction  peak  appears  near  35X  semispan  location 
In  the  Euler-Clebsch  solution;  this  may  be  due  to  the  slope  discontinuity  at 
the  leading-edge  break.  Figure  22  presents  the  solution  at  a  higher  Mach  num- 
ber,  *  0.900  and  a  lower  angle  of  attack,  a  *  2*.  The  Euler-Clebsch  method 
predicts  a  shock  very  close  to  the  trailing  edge,  while  the  ful ly-conservatl ve 
full-potential  method  falls  to  give  a  converged  solution.  Figure  23  presents 
three  solutions  obtained  for  the  F-14  wing/body  using  the  Euler-Clebsch  method 
at  Mach  numbers  M  >  0.85,  0.90  and  0.95,  respectively,  and  at  angle  of  attack 
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Figure  18.  A  cross-sectional  grid  distribution  fo'  the  F-14  wing  body 


a  ^  2°  .  Strong  shocks  are  predicted  on  the  upper  surface  at  all  freestream 
Hach  numbers,' and  a  second  shock  Is  developed  on  the  lower  surface  for  the  M  = 

flD 

0.95  case.  Because  of  the  appearance  of  the  lower-surface  shock  and  a  signif¬ 
icantly  reduced  pressure  plateau  on  the  upper  surface,  the  total  lift  drops  as 
the  freestream  Mach  number  Increases  from  0.900  to  0.950.  The  Kutta  condition 
for  the  solution  of  M^  .  0.95  near  the  wing  .-tip  Is  not  satisfied  exactly  due  to 
the  limitation  of  the  approximate  Euler-Clebsch  assumptions.  As  explained  In 
Section  5.4,  the  Kutta  condition  li  satisfied  under  the  assumption  that  the 
velocity  Jump  across  the  wake,  flq.  Is  much  smaller  than  the  freestream  vel¬ 
ocity.  As  the  shock  becomes  stronger,  the  above  assumption  1s  no  longer  valid; 
therefore  the  pressures  at  the  upper  and  lower  trailing  edges  become  further 
apart,  as  shown  In  the  solutions  of  M  -  0.95.  Whether  the  Kutta  condition 
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Figure  19.  Comparison  of  F-14  wing  alone  and  wing/body  solutions  at  K«,  =  0.800  and 
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can  be  better  satisfied  by  ysinq  the  exact  Euler-Clebsch  method  in  the  three- 
dimensional  case  should  be  further  Investigated. 

The  F-14  wing/body  solutions  presented  in  Figs.  19  through  23  are  obtained  for 
a  wing  with  a  leading-edge  break,  located  at  about  34X  semispan  location,  where 
the  leading-edge  sweep  angle  changes  from  67*  to  22*.  The  tbtal  lift  coef¬ 
ficients  computed  for  this  configuration  at  M  »  0.8  and  a  »  -4.3®,  -2.65®, 

flD 

-1®,  0®,  1®,  2®,  3.05®  and  4.05®  are  presented  In  Table  l  and  Fig.  24  and  com¬ 
pared  with  the  test  data  of  Bavltz  [34].  In  the  calculations,  the  fuselage  is 
extended  to  downstream  infinity  with  constant  cross-section,  as  shown  In  Fig. 
17,  the  tall  section  Is  not  modeled,  and  the  surface  pressure  Integration  for 
the  total  lift  Is  performed  on  the  entire  wing  and  part  of  the  fuse'age  from 
the  nose  to  about  one  root  chord  length  downstream  of  the  wing  trailing  edge. 
The  reference  wing  area  is  chosen  to  be  the  extended  wing  planform  area  which 
Is  coriputed  by  linearly  extending  the  22®  leading-edge  line  and  the  tralllng- 
edge  line  from  the  leading-edge  break  to  the  vertical  symmetry  plane.  The 
computed  total  lift  coefficients  agree  fairly  well  with  the  test  data,  as 
shown  In  Fig.  24,  despite  the  tall  section  and  the  viscous  effect  not  being 
Included  in  the  calculation. 


'ab'e  1.  Comparison  of  experimental  and  computed  total  lift 
coefficients  for  the  F-14  wing/body  conf Iguratlon 
at  M.  --  0.800 


Data  Points 

a 

L  exp 

L'calc 

i 

-4.3® 

-0.410 

-0.3685 

2 

-2.65® 

-0.185 

-0.1 755 

3 

-1* 

0.020 

0.0206 

4 

0® 

0.150 

0.1412 

5 

1® 

'  0.280 

0.2631 

6 

2® 

0.410 

0.3863 

7 

3.05® 

0.515 

0.5183 

8 

4.05' 

0.625 

0.6390 
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7.0  CONCLUDING  REMARKS 


The  ful 1 -potential  and  time-marching  Euler  methods  are  two  of  the  most  promis¬ 
ing  transonic  computational  methods  at  the  present  time.  Each  method  has  Its 
own  advantages  over  the  other  and,  therefore,  has  Its  own  special  applications. 
The  work  reported  has  demonstrated  that  the  full -potential  method  for  transonic 
flowfleld  calculations  can  be  further  Improved  by  Including  the  rotational 
effect.  Introduced  due  to  entropy  Jump  across  the  shocks,  as  In  the  present 
two-  and  three-dimensional  Euler  correction  methods  based  on  the  Clebsch  trans¬ 
formation.  The  approximate  Euler-Clebsch  method  has  been  developed  for  both 
two-  and  three-dimensional  flows,  while  the  exact  Euler-Clebsch  method  Is 
Implemented  only  for  the  tw'‘-d1mens1onal  case.  The  results  for  transonic  air¬ 
foil  flows  show  that  the  approximate  and  the  exact  Euler-Clebsch  solutions 
generally  agree  well.  The  total  pressure  loss  across  the  shocks  Is  Included 
In  the  present  method  1n  conjunction  with  the  use  of  the  Rankine-Hugonlot 
relation,  while  this  Is  neglected  In  the  conventional  full-potential  methods. 
The  convection  of  the  vortldty  Is  computed  along  streamlines  In  the  exact 
Euler-Clebsch  method  and  along  the  mesh  lines  In  the  approximate  Euler-Clebsch 
method,  while  the  vortldty  Is  dissipated  In  the  time-marching  Euler  method 
because  of  the  need  to  add  artificial  dissipation  terms.  The  ‘mplementation 
of  the  present  method  Is  straightforward.  The  present  scheme  based  on  the 
Clebsch  transformation  can  be  applied  In  most  of  the  existing  full-potential 
methods,  not  only  for  calculat'ng  more  accurate  shocxs,  but  also  for  modeling 
two-  and  three-dimensional  rotational  flows  such  as  flows  downstream  of  actu¬ 
ator  disks  [35,36]. 

A  study  of  finite-volume  and  finite-difference  formulations  of  the  full- 
potential  equation  Is  presented.  The  solutions  obtained  by  using  both  formu¬ 
lations  agree  well;  the  finite  volume  method  ^calculates  the  density  term 
Implicitly  while  the  finite-difference  method  does  not.  The  present  method 
computes  transonic  flows  more  accurately  than  the  full-potential  methods  and 
gives  solutions  agreeing  better  with  the  time-marching  Euler  solutions.  How¬ 
ever,  for  flows  with  very  strong  shocks,  the  present  approximate  Euler-Clebsch 
method  does  not  exactly  satisfy  the  Kutta  condition.  Further  work  to  Improve 
the  two-dimensional  exact  Euler-Clebsch  method  and  to  extend  It  to  three 
dimensions  Is  recommended. 
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